Modeling Smectic Layers in Confined Geometries: Order Parameter and Defects 
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We identify a problem with the standard complex order parameter formalism for smectic-A (SmA) 
liquid crystals, which cannot describe certain types of physical defects, and discuss possible alterna- 
tive descriptions of smectic order. In particular, we suggest an approach based on the real smectic 
density variation rather than a complex order parameter. This approach gives reasonable numerical 
results for the smectic layer configuration and director field in sample geometries, and can be used 
to model smectic liquid crystals under nanoscale confinement for technological applications. 



For over forty years, theoretical understanding of smec- 
tic liquid crystals has been based on the complex order 
parameter '^(r) introduced by de Gennes [l] [2], which 
represents the magnitude and phase of layer ordering. 
During this time, the order parameter has been useful in 
many ways. It demonstrated an analogy between smec- 
tic liquid crystals and superconductors, allowing meth- 
ods of solid-state physics to be applied to liquid-crystal 
science [HE]. It led to theories for the nematic-SmA 
and isotropic-SmA transitions, which are strongly af- 
fected by nematic order fluctuations [IHH]. It further 
led to prediction of twist-grain-boundary phases, liquid- 
crystal analogues of the Abrikosov flux lattice in type-II 
superconductors [9]. Most recently, it has led to calcula- 
tions for smectic layer configurations in confined geome- 
tries [ToHin], which may be useful for design of smectic 
devices [17]. 

The purpose of this paper is to point out a problem 
with the complex order parameter description, which af- 
fects some but certainly not all of the work that has been 
done with it. The complex order parameter is unable to 
describe certain types of defects that can realistically oc- 
cur and hence, in some versions of the theory, it predicts 
unphysical configurations of smectic layers. We discuss 
theoretical approaches to resolve this problem, and use 
one of these approaches to calculate smectic layer config- 
urations. 

To see the problem, consider a disclination defect as in 
Fig. [ij In this figure, every point on the plane has a local 
density p(r) = po + ^p{'^)^ with bright and dark regions 
corresponding to higher and lower density, respectively. 
To use the complex order parameter '^(r), we must write 
the local density variation, compared with the average po, 
as Sp{t) = i?e[?/^(r)]. However, it is impossible to asso- 
ciate a unique complex number ip with each point around 
the defect. If we try to make this association, then we 
must say that the phase of tjj increases downward in the 
lower-left quadrant, outward in the right half, upward in 
the upper-left quadrant, and eventually we reach an in- 
consistency. There must be a branch cut where ip changes 
to the complex conjugate ?/;*, as illustrated by the dotted 
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FIG. 1. Disclination in a two-dimensional smectic phase. 
Bright and dark regions correspond to higher and lower den- 
sity, respectively, and the red dotted line is the branch cut in 
ip and fi. 



line. This situation is similar to the well-known problem 
of describing nematic order with a unit vector n(r): go- 
ing around a half-charge disclination, there must be a 
branch cut where n changes to — n. 

The branch cuts in tp and fi occur for the same phys- 
ical reason: neither of these quantities gives an exactly 
correct description of the symmetry of the phase. The 
nematic phase has orientational order along the axis rep- 
resented by ±n, which can be described correctly by a 
tensor. The vector n is often adequate as an approxi- 
mate description, but the sign of fi does not correspond 
to anything physical. Hence, a branch cut in the sign of 
fi is not a physical defect, and it cannot cost any free 
energy. Likewise, the smectic phase has higher density 
at some positions and lower density at other positions, 
and this density variation can be described correctly by 
the real number Sp or Re(ip). The full complex number 
ip may be mathematically convenient as an approximate 
description, but Im{ilj) does not correspond to anything 
physical. Hence, a branch cut in Imi^p) is not a physical 
defect, and it cannot cost any free energy. 

Does this issue with ip affect any calculations? We 



can see two situations with no problem. The first safe 
situation is if the nematic order is well aligned, so that 
there are no half-charged disclinations. For example, in 
a model for the nematic-smectic transition, the nematic 
order might be strong on both sides of the transition, 
with only small spin- wave fluctuations and no topological 
defects. In that case, i/j can describe the smectic order. 

The second safe situation is if branch cuts in tp and 
n compensate for each other in the free energy. For ex- 
ample, de Gennes' original free energy density [1 can be 
written as 

/ = r|^|2 + ||^|^ + C|(V - i</n)^p + fN, (1) 

where q is the favored wavevector of smectic order and 
/at the nematic free energy density. This expression has 
a symmetry on changing (V^, n) -^ (V^*, — n), although it 
does not have symmetries on ip ^ t/j^ or n ^ — n sepa- 
rately. Going around the disclination of Fig. [l] both ip 
and n have branch cuts, where n is discontinuous and 
tp has a discontinuous derivative. However, if the same 
branch cut is used for both variables, then the discon- 
tinuity in h compensates for the discontinuity in \/tp in 
the derivative term of the free energy. As a result, the 
free energy is well-behaved with no singularity across the 
branch cut. 

By contrast, in more recent theories of the isotropic- 
smectic transition [BHS] , n is replaced by the tensor order 
parameter Qij = S{^ninj — \5ij)^ and the free energy 
density is written as 
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In principle, working with Qij instead of n should be 
an improvement, because Qij is uniquely defined every- 
where. However, in this type of theory, there is no way 
to compensate for the branch cut in ip. Hence, the free 
energy is singular along the branch cut, and it cannot 
describe physical defects as in Fig. [l] As a result, there 
is a risk that it might predict unphysical defects instead. 
We can suggest three possible approaches to resolve 
this problem: 

(1) We can go back to de Gennes' original free energy 
in terms of n instead of Qij. In this approach, we must 
make sure that all physical quantities respect the symme- 
try {ip^rv) -^ (V^*,— n), and have no singularity along a 
branch cut where tp and n both switch. For example, in a 
lattice or finite-element version of the model, we must be 
able to switch {ip, n) -^ (^Z^*, — n) on any single site, with- 
out making the switch on neighboring sites. In principle 
this calculation should be possible, but it is awkward to 
work with two unphysical quantities that compensate for 
each other, and hence we will not pursue it further. 

(2) We can replace n by Qij, and replace '^ by a new 
type of mathematical object with the right symmetry. 



This is essentially the approach used by Alexander et 
al. [IHl |19] , who represent smectic order by a height func- 
tion which includes both increasing and decreasing com- 
ponents simultaneously. Alternatively, we might repre- 
sent smectic order by '^ = Re{ip) + i\Im{ip)\^ which can 
be uniquely defined everywhere, with no branch cut. This 
calculation should be possible, but it is awkward to de- 
fine a new type of mathematical object and work out the 
appropriate calculus for it. Hence, in this study we will 
not pursue that approach either. 

(3) We can replace n by Qij, and replace ip by the 
physical density variation Sp. In that case, the theory 
becomes a form of density functional theory, analogous 
to early work on smectic phases [20', l2T] . We find this 
approach the most physical and will use it for the rest of 
this article. 

For this theory, we use a free energy density for the 
SmA phase of the form 
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where dpir) is the local deviation from the average den- 
sity po- This model has a transition from the nematic 
phase when a is above a threshold to the SmA phase 
when a is below the threshold. In the SmA phase, the 
free energy minimum is approximately a sinusoidal den- 
sity wave with wavelength 27: /q. This result is consistent 
with the de Gennes theory. 

We should make four remarks about this free energy. 
First, n only enters through the second-rank tensor niTij^ 
which corresponds to Qij in Refs. [BHH]- Hence, the free 
energy depends only on niUj and 5p^ which are both 
physical, single- valued functions, with no need for branch 
cuts. Second, the free energy includes a third-order term 
of 3p^ . This term is allowed because there is no symmetry 
between high density {5p > 0) and low density {Sp < 0). 
Third, it could include other terms permitted by symme- 
try, such as |V^/)p and (fi • \/Sp)'^. These terms shift the 
nematic-SmA transition and the wavelength of the smec- 
tic density modulation, but do not change the general 
physics discussed here, so we will not consider them fur- 
ther. Fourth, it includes the nematic free energy density 
/at. We use the simplest approximation /at = \K{dinjY 
with a single Frank elastic constant, although it could be 
generalized to different Frank constants. 

In the rest of this paper, we use the free energy (|3| 
to calculate the director orientation and layer configura- 
tion in two-dimensional SmA liquid crystals. We con- 
sider a SmA phase confined inside a small domain with 
boundary conditions requiring defects, and determine the 
lowest-energy structure. These calculations demonstrate 
the feasibility of our density-functional theory, and may 
be useful for designing devices with smectic liquid crys- 
tals under nanoscale confinement [17|. 



To calculate the ground state, we numerically minimize 
the free energy using Monte Carlo simulated annealing. 
The system is discretized using a quasi-regular uniform 
Delaunay triangulation. Its state is described by the di- 
rector orientation and density at each vertex. The total 
free energy is calculated as a sum of free energies of all 
Delaunay cells. The free energy of each cell depends on 
first partial derivatives of the director and second par- 
tial derivatives of density with respect to position. The 
required derivatives are approximated by linear combi- 
nations of the variables at neighboring vertices. 

To calculate the director and its first derivatives in each 
cell, we linearly interpolate each director component be- 
tween the three vertices. The director at the cell's cen- 
ter is determined by interpolation and then normalized 
to unit length. To calculate the density and its second 
derivatives in each cell, we assume the density depends 
quadratically on the coordinates. Interpolation of the 
quadratic function between the closest six vertices gives 
the required derivatives and the density at the center. 

In the Monte Carlo simulated annealing, trial director 
and density changes at any vertex change the free energy 
in all neighboring Delaunay cells, and in cells that use 
the vertex for calculating density derivatives. Thus, for 
each change of director or density at a vertex, the free 
energy is recalculated for several neighboring cells. 

To analyze the numerical results, we want to describe 
the degree of smectic order as a function of position. For 
this purpose, we need the magnitude of density modula- 
tion at the wavelength corresponding to smectic layers. 
The simplest representation for this order parameter is a 
local Fourier transform of the density near any point, at 
the wave vector q along the local director. 
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This quantity is calculated by numerical integration and 
presented as the smectic order parameter in the figures 
below. 

For an initial simulation to illustrate the difference be- 
tween working with Sp{r) and working with '^(r), we 
consider a geometry with two half-charged disclinations. 
In this initial simulation, we assume the director field is 
held fixed and calculate the resulting smectic layer con- 
figuration. When we do this calculation using the free 
energy ([s]), expressed in terms of Sp{r)^ we obtain the 
structure shown in Fig. ^p. This structure is consistent 
with all the symmetries of the SmA phase. The layers are 
equally spaced and normal to the director everywhere. 
The region between the disclinations is a well-ordered 
smectic phase with no line defect. 

By comparison, to illustrate the potential problems 
with using '^(r) and (5ij(r), we perform analogous simu- 
lations using the free energy density from Ref. [6 , similar 




FIG. 2. Simulation of smectic layers using a fixed director field 
(shown by short lines) with two half-charged disclinations. 
(a) Results using the free energy ([5|, expressed in terms of 
'0(r), with parameters a = —1, /3 = 100, bi = —3, 62 = 5, 
ei = — 61 — 852(7r/10)^. (b) Results using the free energy (|3|), 
expressed in terms of Sp{r), with parameters a = —0.1, 6 = 0, 
c = 10, ^ = 27r/10, B = 0.1/^^. 
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As a simplification, we assume here that the magnitude 
of nematic order is constant, so that Qij = fn^nj — |^ij, 
and use the single elastic constant approximation. We 
perform these simulations on a square lattice with deriva- 
tives approximated by standard finite differences. The 
results are shown in Fig. [2^. Here, the color indicates 
\iIj\ (with purple and blue representing higher and lower 
order, respectively), while the brightness indicates Re{ilj) 
(giving the local density). Note that this simulation 
shows a line defect connecting the two disclinations. This 
line defect is a specific example of the branch cut dis- 
cussed earlier: a boundary where V?/^ -^ Vifj* . As an 
artifact of the model, this boundary has a free energy 
penalty, which is linearly proportional to the distance be- 
tween disclinations and hence binds the disclinations to- 
gether. To minimize the free energy, the algorithm shifts 
some of the penalty from /m('0) to i?e('^), and hence it 
appears as a line defect in the density wave. That line 
defect is unphysical, because there is no need for any line 
defect there; the density wave can have a smooth curva- 
ture as calculated through density functional theory. 

For a second example, we perform simulations where 
the director field and layer configuration can both relax. 
We consider the circular geometry shown in Fig. [3J with 
boundary conditions on the director requiring a single 
half-charged disclination. Subject to that constraint, the 
director and layers relax together inside the domain. 

When we do the calculation with the free energy p|), 
using ^p(r), we obtain the structure shown in Fig. [3p, 
which has a single half-charged disclination at the center. 
The smectic layers form a relaxed configuration about 
the disclination. There is a point defect in the layers at 




FIG. 3. Simulations of a circular domain with boundary con- 
ditions requiring a single half- charged disclination. (a) Re- 
sults using the free energy ([5]), with parameters a = —1, 
13 = 100, bi = -3, 62 = 5, ei = -bi - 862(7^/10)^ 
K = 0.0025. (b) Results using the free energy (Is]), with pa- 
rameters a = -0.1, 6 = 0, c = 10, ^ = 27r/10, B = 0.1/^^, 
K = 0.008. 



the disclination core, where we can see a reduction in 
the smectic order parameter defined by the local Fourier 
transform. Everywhere else, the layers are well-ordered 
and equally spaced. 

By comparison, when we do the calculation with the 
free energy ([5|, using V^(r), we obtain the structure shown 
in Fig. [3^. Once again, we see a line defect in the layers 
coming out of the disclination. Because there is no other 
disclination where the line defect can terminate, it runs 
all the way to the boundary. This line defect is not re- 
quired by the symmetry of the smectic phase; it is just 
an artifact of the complex order parameter formalism. 

For further examples of the density functional theory, 
we perform simulations of the circular domains shown in 
Fig. [4J Here, the director field has tangential boundary 
conditions (parts a and b) or radial boundary conditions 
(c and d). In either case, it must have total topological 
charge of +1. The density modulation has free boundary 
conditions. We use two values of the Frank elastic con- 
stant K compared to the nematic-smectic coupling 5, 
and hence two values of the length scale A = (K/B)^^'^. 
This characteristic smectic length scale is large in a and 
c, and smaller in b and d. In all cases, free energy min- 
imization gives a configuration with two disclinations of 
topological charge +1/2 each, not a single disclination 
of +1. In the two cases with high A, the director has a 
smooth variation between the disclinations, and the lay- 
ers adapt to the director, with small variations in the 
layer spacing. In the two cases with smaller A, the lay- 
ers are equally spaced over most of the domain, and 
the director adapts to the layers, with director varia- 
tion concentrated in small regions near the boundary. 
These results are physically reasonable, and correspond 
to what might be observed for smectic liquid crystals un- 
der nanoscale confinement [17]. 

As a final point, we note that the complex order param- 
eter formalism has a second problem, albeit less severe 




FIG. 4. Simulations of a circular domain with boundary con- 
ditions requiring tangential alignment (parts a and b) or ra- 
dial alignment (c and d) of the director. The density func- 
tional free energy Q is used. In a and c, parameters are 
a = -5, 6 = 0, c = 5, B = 10"^ q = 40, K = 0.3. In b and 
d, the Frank constant is reduced to K = 0.05. 



than the line defects discussed above. In this formalism 
all positions are equivalent, with the same free energy 
density, regardless of whether they are density maxima, 
minima, or anywhere in between. For that reason, a dis- 
location in the smectic layers is equally likely to occur 
anywhere within the layer structure. As a result, there 
is no Peierls-Nabarro energy barrier for dislocation glide 
from layer to layer, the process illustrated in Fig. 9.17 
of Ref. [22J. Indeed, if we multiply '^(r) by a uniform 
phase factor e*"^, we explicitly shift the layers with re- 
spect to the defect position, and obtain dislocation glide 
at zero energy cost. This symmetry is unphysical: in the 
real system, the density maxima and minima are special 
positions, and there must be some free energy difference 
between defects at those positions and defects at other 
positions. The density formalism shows these special po- 
sitions explicitly, and hence predicts a nonzero Peierls- 
Nabarro barrier for dislocation glide [23 . 

In conclusion, we have identified a problem with the 
complex order parameter formalism for modeling smec- 
tic liquid crystals, which affects calculations that should 
be able to produce half-charged disclinations. We have 
discussed three ways to solve this problem: canceling the 
branch cut in the smectic order parameter with another 
branch cut in the nematic director, defining a new type of 
mathematical object for the smectic order parameter, or 
working with the physical density variation Sp. All three 



of these approaches should be possible, but the third ap- 
proach seems to us the simplest and most physical, and it 
has the further advantage of a nonzero Peierls-Nabarro 
barrier. We have demonstrated that it works for sam- 
ple geometries, and it has potential for future design of 
nanoscale smectic devices. 
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SUPPLEMENTAL MATERIAL 



This note presents explicit numerical calculations of 
the energy of an edge dislocation in the 2D SmA phase, 
using the two theoretical approaches discussed in the Let- 
ter. In the complex order parameter formalism, the dis- 
location energy is independent of position with respect to 
the layer structure, and hence there is no Peierls-Nabarro 
energy barrier for dislocation glide. By contrast, in the 
real density formalism, the dislocation energy depends 
on position and the Peierls-Nabarro barrier is nonzero, 
as is physically realistic. 

In the 2D (x, y) plane, the local density is given by 
p{x,y) = po -\- Sp{x^y)^ which is related to the complex 
order parameter by Sp{x^y) = Re[ilj{x,y)]. To describe 
a single edge dislocation, we choose the complex order 
parameter ^l){x^y) = e*('?^+*^(^'^)+^^), where ^{x^y) = 
arg(x + iy) is the phase variation and A^ is a constant 
phase offset. The director n{x^y) is chosen as a unit 
vector along the gradient of '^(x, ?/). 

Figure 1 shows a visualization of ^{x^y); the branch 
cut starts at (x = 0, ?/ = 0) and goes in positive x di- 
rection. The first column of pictures in Fig. 2 shows cor- 
responding visualizations of the density variation around 
the dislocation. Note that the constant phase offset A$ 
defines the position of the dislocation with respect to the 
layer structure. When A<I> = the dislocation occurs at 
a density minimum (shown in blue); when A<I> = tt it 
occurs at a density maximum (shown in red). For inter- 
mediate A$ it occurs at a lower-symmetry point between 
those extremes. 

In the complex order parameter formalism, the free 
energy density is given by 



fw = rm 



' + ^\i>\' + C\{V-iqh)^\'. 



(1) 




FIG. 1. Visualization of the phase variation ^(x,y) around a 
single edge dislocation. 



As a specific example, the second column of pictures in 
Fig. 2 shows the free energy density for the parameters 
r = —5, u = b, q = 1^ and C = 1/q^. From these 
pictures, we can see that the free energy density is sharply 
peaked at the dislocation, and decays rapidly away from 
the dislocation. The free energy density of Eq. (1) is 
clearly independent of the constant phase offset A$, and 
indeed all the pictures in this column are identical. As a 
result, the integrated free energy is also independent of 
A<I>. This result implies that the dislocation can move 
with respect to the layer structure, from row to row in 
the table, with no energy cost. In the terminology of 
dislocation theory, we would say that the Peierls-Nabarro 
energy barrier for dislocation glide is zero. This result is 
physically unrealistic, and shows that the model has more 
symmetry than the actual SmA phase. 

By comparison, in the real density formalism, the free 
energy density is given by 

f{5p) = \bp^ + \bp^ + I J/ + B [(5,5,- + q^n,ni)bp\\ 

(2) 
For specific examples, the third and fourth columns of 
pictures in Fig. 2 show the free energy density calculated 
for the parameters a = — 10, c = 10, g' = 1, 5 = 0.1/g^, 
and 6 = (in the third column) and 6 = 1 (in the fourth 
column). The h term is important because it is the only 
term considered here that is odd in (5p, and hence the 
only term that distinguishes between density minima and 
maxima. From these images, we can make several obser- 
vations. First, the free energy density is not uniform 
but periodic in the smectic layer structure. If 6 = 0, 
there are equal free energy valleys at the density minima 
and maxima. If 6 > 0, the symmetry between minima 
and maxima is broken (as is physically realistic), and the 
deepest free energy valleys are at the density minima. 
Furthermore, there is additional free energy associated 
with the dislocation itself. Most importantly, the free 
energy plots change as the constant phase offset A<I> is 
varied, i. e. as the dislocation moves with respect to the 
layer structure. Hence, this model does not have the un- 
physical symmetry found in the complex order parameter 
formalism. 

To calculate the barrier for dislocation motion, we 
must integrate the free energy density to find the to- 
tal free energy as a function of A<I>. In this calcula- 
tion, it is important to integrate over an integer num- 
ber of layers, so that the result is not influenced by the 
number of fractional layers within the integration region. 
Hence, we define the integration region by \qx\ < IOtt 
and \qy + <I>(x, y)\ < Stt, as shown in Fig. 3. 

Figure 4 presents graphs of the total free energy as a 
function of A<l>. We can see that it has a periodic series of 
peaks and valleys as the dislocation moves with respect to 
the layer structure. If 6 = 0, the valleys occur whenever 
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FIG. 2. Behavior as the dislocation is displaced with respect to the layer structure. Plots show the smectic density variation, 
the free energy density for the complex order parameter formalism, and the free energy density for the real density formalism 
(calculated for the cubic parameter 6 = and 6/0). Red represents highest values of density or free energy; blue represents 
lowest values. 
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FIG. 3. Integration of the free energy density over a region containing an integer number of smectic layers. 
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FIG. 4. Integrated free energy in the real density formalism as the dislocation is displaced with respect to the layer structure, 
as given by the parameter A$. (a) Parameters a = —10, c = 10, ^ = 1, B = 0.1/^^, and 6 = 0, so that there is a symmetry 
between density minima and maxima, (b) Parameter 6=1, breaking the symmetry between density minima and maxima. 



A<^ is a multiple of tt, i. e. whenever the dislocation valleys at density maxima are less deep. The Peierls- 

is at either a density minimum or maximum. If 6 > 0, Nabarro energy barrier for dislocation glide is then the 

the deepest valleys occur when A$ is a multiple of 27r, difference in free energy between the deepest valleys and 

i. e. when the dislocation is at a density minimum; the highest peaks in this plot. 



